module testmomsrot
implicit none
contains

function eubeta3(beta,eta,si)
!E(epstar^3)
use ghs
double precision beta,eta,si,eubeta3,nu,gam,c
double precision h1,h2,h3,p1,p2

nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
h1=1.0d0
h2=ehr(nu,gam,2)
h3=ehr(nu,gam,3)

p1=h3-1.0d0-3.0d0*h2+3.0d0*h1 !E((h-1)^3)
p2=h2-h1 !E(h*(h-1))

eubeta3=(c*beta)**3.0d0*p1+3.0d0*c**2.0d0*beta*p2
end function eubeta3

function eubeta4(beta,eta,si)
!E(epstar^4)
use ghs
double precision beta,eta,si,eubeta4,nu,gam,c
double precision h1,h2,h3,h4,p1,p2,p3

nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
h1=1.0d0
h2=ehr(nu,gam,2)
h3=ehr(nu,gam,3)
h4=ehr(nu,gam,4)

p1=h4-4.0d0*h3+6.0d0*h2-4.0d0*h1+1.0d0 !E((h-1)^4)
p2=h3-2.0d0*h2+h1 !E(h*(h-1)^2)
p3=h2
eubeta4=(c*beta)**4.0d0*p1+6.0d0*c**3.0d0*beta**2.0d0*p2+3.0d0*c**2.0d0*p3
end function eubeta4

function eubeta5(beta,eta,si)
!E(epstar^5)
use ghs
double precision beta,eta,si,eubeta5,nu,gam,c
double precision h1,h2,h3,h4,h5,p1,p2,p3,p4

nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
h1=1.0d0
h2=ehr(nu,gam,2)
h3=ehr(nu,gam,3)
h4=ehr(nu,gam,4)
h5=ehr(nu,gam,5)

p1=h5-5.0d0*h4+10.0d0*h3-10.0d0*h2+5.0d0*h1-1.0d0 !E((h-1)^5)
p2=h4-3*h3+3*h2-h1 !E(h*(h-1)^3)
p3=h3-h2  !E(h^2*(h-1))
eubeta5=(c*beta)**5.0d0*p1+10.0d0*c**4.0d0*beta**3.0d0*p2+15.0d0*c**3.0d0*beta*p3
end function eubeta5

function eubeta6(beta,eta,si)
!E(epstar^6)
use ghs
double precision beta,eta,si,eubeta6,nu,gam,c
double precision h1,h2,h3,h4,h5,h6,p1,p2,p3,p4,p5

nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
h1=1.0d0
h2=ehr(nu,gam,2)
h3=ehr(nu,gam,3)
h4=ehr(nu,gam,4)
h5=ehr(nu,gam,5)
h6=ehr(nu,gam,6)


p1=h6-6.0d0*h5+15.0d0*h4-20.0d0*h3+15.0d0*h2-6.0d0*h1+1.0d0 !E((h-1)^6)
p2=h5-4.0d0*h4+6.0d0*h3-4.0d0*h2+h1 !E(h*(h-1)^4)
p3=h4-2.0d0*h3+h2 !E(h^2*(h-1)^2)
p4=h3 !E(h^3)
eubeta6=(c*beta)**6.0d0*p1+15.0d0*c**5.0d0*beta**4.0d0*p2+45.0d0*c**4.0d0*beta**2.0d0*p3&
&+15.0d0*c**3.0d0*p4
end function eubeta6

function eubeta7(beta,eta,si)
!E(epstar^7)
use ghs
double precision beta,eta,si,eubeta7,nu,gam,c
double precision h1,h2,h3,h4,h5,h6,h7,p1,p2,p3,p4,p5

nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
h1=1.0d0
h2=ehr(nu,gam,2)
h3=ehr(nu,gam,3)
h4=ehr(nu,gam,4)
h5=ehr(nu,gam,5)
h6=ehr(nu,gam,6)
h7=ehr(nu,gam,7)

p1=h7-7.0d0*h6+21.0d0*h5-35.0d0*h4+35.0d0*h3-21.0d0*h2+7.0d0*h1-1.0d0 !E((h-1)^7)
p2=h6-5.0d0*h5+10.0d0*h4-10.0d0*h3+5.0d0*h2-h1 !E(h*(h-1)^5)
p3=h5-3.0d0*h4+3.0d0*h3-h2 !E(h^2*(h-1)^3)
p4=h4-h3 !E(h^3*(h-1))
eubeta7=(c*beta)**7.0d0*p1+21.0d0*c**6.0d0*beta**5.0d0*p2+105.0d0*c**5.0d0*beta**3.0d0*p3&
&+105*c**4.0d0*beta*p4
end function eubeta7

function eubeta8(beta,eta,si)
!E(epstar^8)
use ghs
double precision beta,eta,si,eubeta8,nu,gam,c
double precision h1,h2,h3,h4,h5,h6,h7,h8,p1,p2,p3,p4,p5

nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
h1=1.0d0
h2=ehr(nu,gam,2)
h3=ehr(nu,gam,3)
h4=ehr(nu,gam,4)
h5=ehr(nu,gam,5)
h6=ehr(nu,gam,6)
h7=ehr(nu,gam,7)
h8=ehr(nu,gam,8)

p1=h8-8.0d0*h7+28.0d0*h6-56.0d0*h5+70.0d0*h4-56.0d0*h3+28.0d0*h2-8.0d0*h1+1.0d0 !E((h-1)^8)
p2=h7-6.0d0*h6+15.0d0*h5-20.0d0*h4+15.0d0*h3-6.0d0*h2+h1 !E(h*(h-1)^6)
p3=h6-4.0d0*h5+6.0d0*h4-4.0d0*h3+h2 !E(h^2*(h-1)^4)
p4=h5-2.0d0*h4+h3 !E(h^3*(h-1)^2)

eubeta8=(c*beta)**8.0d0*p1+28.0d0*c**7.0d0*beta**6.0d0*p2+210.0d0*c**6.0d0*beta**4.0d0*p3&
&+420.0d0*c**5.0d0*beta**2.0d0*p4+105.0d0*c**4.0d0*h4
end function eubeta8

!!!
function eub1hn(beta,eta,si,n)
!E(ubeta^1*h^n)
use ghs
integer n
double precision beta,eta,si,eub1hn,nu,gam,c
nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
eub1hn=c*beta*(ehr(nu,gam,n+1)-ehr(nu,gam,n))
end function eub1hn


function eub2hn(beta,eta,si,n)
!E(ubeta^2*h^n)
use ghs
integer n
double precision beta,eta,si,eub2hn,nu,gam,c
double precision h0,h1,h2,p1,p2
nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
h0=ehr(nu,gam,n)
h1=ehr(nu,gam,1+n)
h2=ehr(nu,gam,2+n)

p1=h2-2.0d0*h1+h0 !E(h^n(h-1)^2)
p2=h1 !E(h^(n+1))

eub2hn=(c*beta)**2.0d0*p1+c*p2
end function eub2hn

function eub3hn(beta,eta,si,n)
!E(ubeta^3*h^n)
use ghs
integer n
double precision beta,eta,si,eub3hn,nu,gam,c
double precision h1,h2,h3,h0,p1,p2
nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
h0=ehr(nu,gam,n)
h1=ehr(nu,gam,1+n)
h2=ehr(nu,gam,2+n)
h3=ehr(nu,gam,3+n)

p1=h3-h0-3.0d0*h2+3.0d0*h1 !E(h^n(h-1)^3)
p2=h2-h1 !E(h^(n+1)*(h-1))

eub3hn=(c*beta)**3.0d0*p1+3.0d0*c**2.0d0*beta*p2
end function eub3hn

function eub4hn(beta,eta,si,n)
!E(ubeta^4*h^n)
use ghs
integer n
double precision beta,eta,si,eub4hn,nu,gam,c
double precision h0,h1,h2,h3,h4,p1,p2,p3
nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
h0=ehr(nu,gam,n)
h1=ehr(nu,gam,1+n)
h2=ehr(nu,gam,2+n)
h3=ehr(nu,gam,3+n)
h4=ehr(nu,gam,4+n)

p1=h4-4.0d0*h3+6.0d0*h2-4.0d0*h1+h0 !E(h^n(h-1)^4)
p2=h3-2.0d0*h2+h1 !E(h^(n+1)*(h-1)^2)
p3=h2 !E(h^(2+n))
eub4hn=(c*beta)**4.0d0*p1+6.0d0*c**3.0d0*beta**2.0d0*p2+3.0d0*c**2.0d0*p3
end function eub4hn

function eub5hn(beta,eta,si,n)
!E(ubeta^5*h^n)
use ghs
integer n
double precision beta,eta,si,eub5hn,nu,gam,c
double precision h0,h1,h2,h3,h4,h5,p1,p2,p3,p4

nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
h0=ehr(nu,gam,n)
h1=ehr(nu,gam,n+1)
h2=ehr(nu,gam,n+2)
h3=ehr(nu,gam,n+3)
h4=ehr(nu,gam,n+4)
h5=ehr(nu,gam,n+5)

p1=h5-5.0d0*h4+10.0d0*h3-10.0d0*h2+5.0d0*h1-h0 !E(h^n(h-1)^5)
p2=h4-3*h3+3*h2-h1 !E(h^(n+1)*(h-1)^3)
p3=h3-h2  !E(h^2*(h-1))
eub5hn=(c*beta)**5.0d0*p1+10.0d0*c**4.0d0*beta**3.0d0*p2+15.0d0*c**3.0d0*beta*p3
end function eub5hn

function eub6hn(beta,eta,si,n)
!E(ubeta^6*h^n)
use ghs
integer n
double precision beta,eta,si,eub6hn,nu,gam,c
double precision h0,h1,h2,h3,h4,h5,h6,p1,p2,p3,p4,p5

nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si
c=ghcfun1(beta*beta,eta,si)
!!
h0=ehr(nu,gam,n)
h1=ehr(nu,gam,n+1)
h2=ehr(nu,gam,n+2)
h3=ehr(nu,gam,n+3)
h4=ehr(nu,gam,n+4)
h5=ehr(nu,gam,n+5)
h6=ehr(nu,gam,n+6)


p1=h6-6.0d0*h5+15.0d0*h4-20.0d0*h3+15.0d0*h2-6.0d0*h1+h0 !E((h-1)^6)
p2=h5-4.0d0*h4+6.0d0*h3-4.0d0*h2+h1 !E(h*(h-1)^4)
p3=h4-2.0d0*h3+h2 !E(h^2*(h-1)^2)
p4=h3 !E(h^3)
eub6hn=(c*beta)**6.0d0*p1+15.0d0*c**5.0d0*beta**4.0d0*p2+45.0d0*c**4.0d0*beta**2.0d0*p3&
&+15.0d0*c**3.0d0*p4
end function eub6hn


!!!!!!!!!!!!!!!
function ehr(nu,gam,k)
!E(ht**k)
use bessels
double precision nu,gam,ehr
integer k,i
ehr=1
do i=1,k
	ehr=ehr*besselr(nu+i-1,gam)/besselr(nu,gam)
end do

end function ehr

end module testmomsrot